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Abstract 

We give a. simple faotorization^of an arbitrary hermitian, positive definite 
matrix in which the factors are well-conditioned, hermitian, and positive 
definite. In fact, given knowledge of the extreme eigenvalues of the original 
matrix A, we can achieve an optimal improvement, making the condition 
numbers of each of the two factors equal to the square root of the condition 
number of A. 

We apply this technique to the solution of hermitian, positive definite 
Toeplitz systems. Large linear systems with hermitian, positive definite 
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Toeplitz matrices arise in some signal processing applications. We give a 
stable fast algorithm for solving these systems that is based on the precondi- 
tioned conjugate gradient method. The algorithm exploits Toeplitz structure 
to reduce the cost of an iteration to O(nlogn.) by applying the fast Fourier 
Transform to compute matrix- vector products. We use our matrix factoriza- 
tion as a preconditioner. 


1 Introduction 

We give a simple factorization of an arbitrary hermitian, positive definite 
matrix A in which the factors are hermitian, positive definite, and are sub- 
stantially better conditioned than the original matrix A. In fact, given knowl- 
edge of the extreme eigenvalues of .4, we can achieve an optimal improve- 
ment, making the condition numbers of each of the two factors equal to the 
square root of the condition number of A. The factorization is of the form 
A — {A + — n(A + /t/) -1 ). We discuss the optimum choice of /i in 

Section 3. 

Consider the linear system Ax = b where 4 is an n x n hermitian, positive 
definite Toeplitz matrix, .4 = [a tJ ] = [« Jt ] - [«| t _y|]. Several direct methods 
for solving such a system using 0(n 2 ) arithmetic operations are known [18], 
[14], [12], [3]. In addition, some newer methods that take 0(nlog 2 n) opera- 
tions have been developed [4], [5], [13], [ll], [2]. The method of Gragg and 
Ammar, for example, requires 8?/ log 2 n real arithmetic operations for a real 
Toeplitz system. Stability of these methods is discussed by Bunch [6]. 

Recently, some new attention has been given to the preconditioned con- 
jugate gradient method as a Toeplitz solver. The motivation is that a single 
iteration of this method can be implemented, using the fast Fourier Trans- 
form, at a cost of 0(n log n) arithmetic operations. The hope is that with 
suitable preconditioners, the number of conjugate gradient iterations can be 
made small enough to make the method practical. One approach has been 
to use circulant approximations to A as preconditioners [16], [7]. These work 
remarkably well for some Toeplitz matrices (finite sections of a singly in- 
finite matrix A^ = [« t _j] for which ££L 0 a* < °°) [8], [17]. As we shall 
show in Section 4, however, they are not always effective. Here we use in- 
stead the general preconditioning strategy discussed above. This strategy is 
most useful when A is Toeplitz, for in this case, the factor, (.4 + fil), is still 
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Toeplitz. Conjugate gradient iterations are therefore inexpensive. Moreover 
its inverse has a representation as the difference between two products of 
triangular Toeplitz matrices (the Gohberg-Semencul formula) or a triangu- 
lar Toeplitz matrix and a circulant matrix (the Ammar-Gader formula); this 
allows iteration with the other factor, (/ — fi(A + pi) *) to be carried out 
cheaply. 

In the next section we give a general outline of our method. In Section 3 
we give the details of this method, providing an optimum shift paiameter g. 
We compare it with several competitors in Section 4. 


1.1 Notation 

For z E C , z" denotes the complex conjugate of z. For a complex matrix 
A , A h denotes the conjugate transpose of A , and A(A) denotes the set of 
eigenvalues of the square matrix A. A matrix A is hermitian, positive definite 
(hereafter h.p.d.) if A = A H and for all a € A(A),a > 0. Let A be li.p.d. 
and let a x be the largest and a„ the smallest of A’s eigenvalues. Then 
k(A) = Qj/a n is called the (spectral) condition number of A. 


2 A condition-improving matrix factoriza- 
tion 

Let A be a given h.p.d. matrix. Let p € R- Let 

B=A + gI (1) 

and 

C = I-yB~ l . (2) 

Lemma 1 Let A be a given h.p.d. matrix. Let B and C be given by (1) and 
(2). Then A = BC = CB. If -p is not an eigenvalue of A then both B and 
C have inverses and A~ l = C~ 1 B~ ] = B l C 

Proof: A = B - gl = B(I - pB~ l ) = {I - pB~ l )B. Invertibility of B is 

obvious. And C = B~ l A must therefore be invertible too. This proves the 
lemma. 
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Let the eigenvalues of A, B , and C be given by 

{^n <•••</?!} = A (B), 
{in < ■ •• < 7l} = A(C). 

By (1) and (2) we have 



+ 

S 

II 

(3) 


7 j — 1 i l Pj • 

(4) 

Lemma 2 Let B andC be 

given by (1) and ( 2 ). 

Then the condition numbers 

of B and C arc given by 

«(B)= a ' + " 
cr„ + n 

(5) 

and 

«i («» + /*) 
k(U) = — 7 r . 

Ctn{ot 1 + H) 

(6) 

3 

S 2 

J/i 

> 

*1 

+*** 

V 

o 

k(A) = k(B)k(C) . 

(?) 


Proof: Immediate. 

Lemma 3 Let n = ^/Q x a n . Then k(B) = k(C) — \J n{A). 
Proof: Immediate. 


In most cases, this factorization does not lead to a reduction in the cost of 
inverting A or of solving Ax = b. For if we employ an algorithm for inverting 
B and C, such as Gaussian elimination, that is insensitive to condition num- 
ber, we might as well have used it on A\ if the cost of this algorithm depends 
on log k, as it does for Newton’s method (see, for example, [15]) then we 
have gained nothing. If the cost exceeds log/c, as it does for most iterative 
methods, we may have gained something. But iteration to solve Cx = z is 
expensive because we need to solve a system with B at each iteration. Some 
form of inner/outer iteration method could be used to reduce this cost. 
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There is a case, however, in which we can solve many systems with B in 
little more time than it takes to solve one. When A is an h.p.d. Toeplitz 
matrix, then so is B. Then we may use the Gohberg-Semcncul formula, or 
a recent variant, the Ammar-Gader formula, to represent B~ l . This reduces 
the cost of one iteration of CG for C to 0(n log n). 


3 A Fast Toeplitz Solver 


Let us apply the matrix factorization of the previous section to the solution 


of a linear system 


Ax = b 


(«) 


where A = is an n x n h.p.d. Toeplitz matrix, a t j — aji — We 

suppose that the extremal eigenvalues of A have been estimated (more on 
this later on, in Section 5). In this case, the matrix B is also h.p.d. and 
Toeplitz, is completely defined by its first column, and its inverse can be 


represented as 


B~ x = L X L* - L 2 L j 2 j 


where L x and L 2 are lower triangular Toeplitz matrices whose elements de- 
pend only on the first column of B~ l ([9], [18]). Recently, an alternative rep- 
resentation by Toeplitz-circulant products was given by Ammar and Gadei 
[ 1 ]: 

B~ x = L x F7 - L 2 E (9) 


where L x and L 2 are again lower triangular Toeplitz, and E is circulant. 
These factors are again functions only of the first column of B~ x . 

This leads to the following algorithm: 


Algorithm 1: 

Input: An h.p.d. Toeplitz matrix A, a vector 6, and a shift parameter \i. 
Output: A -1 b. 

Method: 

(Step 1) Solve By = e u where e t = (1,0, . . . , 0) // . Construct L\,LiiE 
satisfying (9). 

(Step 2) Solve Bz = b. 

(Step 3) Solve Cx = 2 . Return x. 
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We solve the linear systems at Steps 1 and 3 of Algorithm 1 by the 
conjugate gradient (CG) method [10]. At Step 2, we use the representation 
(9) constructed at Step 1; there is no iteration and only FFTs of n-vectors 
are required. 

In the CG method, each iteration costs 5n flops and one matrix-vector 
product. For the Toeplitz matrix B in Steps 1 and 2, we compute the matrix 
vector product Bu by a standard technique of embedding B in a circulant 
matrix of order 2 n and appending n zeros to u; a circulant matrix times a 
vector is a convolution, for which the fast Fourier Transform is used: 

u = F 2b ([«, 0]); 
v — u. * B] 
v = W 2n (v); 

Bu = u(l : n); 

Here 0 represents the zero vector of order ri; Ft and Wk are the forward and 
inverse discrete Fourier Transform operators on C k , the operator .* represents 
elementwise multiplication of two vectors, v(l : n) is the first n components 
of v, and 

B = r. 2n ([b 

Note that £} is real and this reduces the cost of the elementwise multiplica- 
tion. The cost of a CG iteration with B is therefore 

Cost(B) = 4$(n) 

where $(n) is the cost of an n-point FFT. (This does not include the cost of 
computing B, which is computed once and for all.) 

For Step 3 above, we use CG, too. We require the product Cu and hence 
also B~ l u at every iteration. Using the Ammar-Gader formula, this cost may 
be reduced to 

Cost(C) = 7 $(n) 

( see [!])• 

The alternative of using B as a preconditioner in the preconditioned CG 
method is less attractive than Algorithm 1. We would in effect be solving 

Cx = B~ l b 
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by the conjugate gradient method. Thus, we have an equivalent of Steps 
2 and 3 above. But at each iteration, we would have to form Cp for some 
vector p as B~ x Ap. The cost of each iteration would therefore be ll$(n) 
rather than 7 <J>(n) as it is in our implementation. 

3.1 The Optimal Shift 

The choice p = v /a7o^’ is not optimal. Since each iteration of CG with the 
matrix C costs roughly 1.75 times as much as a CG step with B, wc would 
be better off to choose p to make C better conditioned at the expense of 
worsening the condition number of B. It would be reasonable to choose p to 
minimize the estimate of the total work given by 


(4 n B + 7nc)$(n), 

where ng is the number of CG iterations at Step 1 above, nc the number of 
CG iterations at Step 3 above. According to [10], we may model these by 



n B = F\Jk(B) 

(10) 

and 

n C — F yfit{C) 

(11) 

with F a constant. 

Then, by Lemma 2, 



n B n c = F 2 \j n{A) = M = const. 


Therefore, for any constant I< (and motivated by the estimate above we may 
think of I< = 7/4), 


f(n c ) = n B + Kn c 
M , 

= b A n c 

nc 


is minimized at 

[m 

nc = V^’ 

M 

n B = — 

n c 

= Vmk 

= Kn c . 

(12) 
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Now we want to choose the best shift parmeter /x. In view of (10), (11), and 
(12) we choose fi to solve 


k(B) = K 2 k(C), 


(13) 


which becomes, using (5) and (6), 

cr„(a 2 + 2ya l + fj 2 ) = K 2 ax{a 2 n + 2/xa„ + /x 2 ) . 

Now change variables from /x to m, where /x = m^/ot then m satisfies 
m 2 (K 2 ai — a n ) + m(2(K 2 — l)\/ai«T) + (^ 2<1 n — a i) = 0 
or, dividing by a n > 0, with k = n(A) = a ( /or n , 

2 (A' 2 k - 1 ) + m(2(K 2 - l)v^) + (A' 2 - «) = 0 • 


m 


The roots are 


m± = 


~{I< 2 ~1)^±K(k-1) 
IOk - 1 


(14) 


Note that m + > 0 only when k > K 2 . If k < K 2 there is no possibility of 
satisfying (13) anyway, in view of (7). Thus, we assume that k > K 2 and 
that n = m+y/a^ > 0. For « > K 2 the root m + is monotone increasing 
with k and is asymptotic to 1/A'. 


Lemma 4 Let /x be given by my/a. i« n where m — m + (see (14) above). Then 

k(B) = (15) 

k(C) = (16) 

Proof: According to the derivation of (14), the relation (13) holds. Then, 

by (7), 

k = k(B)k(C) 

= (. k(B)/K ) 2 , 

whence (15) follows. And (16) follows from (13) and (15). This proves the 
lemma. 
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Now, assume (10) and (11) and choose y = m +y /a^a^ so that (13) holds. 
In our application, K = 1.75. The total cost is then 


(4n B + 7nc)$(u) = 


4 (n B + Kn c )${n) 
4(2 n B )<fr(n) 

8 y/ K {B)F${n) 
8y/KK l ' A F<i>(n) 
4v/7 k 1/4 F${ti). 


(17) 


For comparison, let n GG be the number of iterations required by CG for /l. 
We assume the model 

Cost(CC) = 4n CG $(n) = 4F$(n)^//c(A) . (18) 


Compared with (18) we can see that the new method is an improvement for 
k > 49. 


3.2 Recursive Preconditioning 

The factorization A = BC can be used recursively. In particular, let us 
consider solving the equation By = Cj at Step 1 of Algor; Jim 1 by using this 
preconditioned conjugate gradient solver. We now solve for two optimum 
shift paramemters; yi, which we use in solving for y in Step 1, is given by 
the theory above in which we substitute /?i , and k(B) for ai,a n , and 
/c(A). The other, /i, which we use to define £?, is now chosen to minimize the 
sum of the cost of Step 1, which by (17) is 

4\/7,c(B) 1/4 F#(») , 

and the cost of Step 3, which is 

7/c(C) 1/2 F$(n) 

subject to (7). The solution is to take 
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and 



Making the substitution // = myJoc\Ot n as before leads to a cubic in rn that 
has a positive root for all k > 49/16. The overall work becomes 

✓ *7 \ 2/3 

s(-J Fi(»)«(A) 1/f . 

This is better that 4F$(n)/c(A) 1/2 (compare (18)) for k(A) > 140 and better 
than 4\/7F$(7i)(«:(/l)) 1/4 (compare (17)) for «:(i4) > 3222. 

4 Numerical Results and Conclusions 

First, we used Algorithm I to solve the system (8) where A was as follows. 
We generated a real time series 

Sk = 2sin (Xk) + sin(4i fc + tt/6) + 

where the n sample points x k were uniformly spaced in [0, 2n) and i/ k was a 
Gaussian random variable (white noise) with mean zero and variance (power) 
a 2 . Next, we took aj to be the autocorrelation of {s/t} with lag j, that is, 

n— 1 

a 3 = S ^ Sk +3 
k—0 

where the index k + j was taken modulo n. 

To choose /*, we used the Lanczos algorithm to approximately tridiago- 
nalize A , then computed the extreme eigenvalues n and r m of the resulting 
m x m symmetric tridiagonal matrix. We used the approximations 


on ~ , 

O' n ~ * 

Tables 1 and 2 give experimental results for a = 10 and <7 = 1. In each 
table we list the number of sample points, n; the number of CG iterations 
at Step 1 of Algorithm 1, ns', the number of CG iterations at Step 3 of 
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n 

n B 

nc 

BHEOI 

nca 

\\Ax - y h 

\\Ax C G-y\\ 

16 

17 

7 

29 

17 

1.8(12) 

1.2(- 13) 

64 

39 

13 

62 

46 

9.7(-12) 

6.9(-12) 

256 

50 

18 

82 

73 

l-2(-ll) 

8. 1 (- 1 2) 

1024 

70 

23 

110 

124 

5.3(- 11) 

2.9(-n) 

4096 

59 

38 

126 

- 

162 

1.2(10) 

6- 1(11) 


Table 1: Results for a = 10 . 


n 

nB 

nc 


nca 

\\Ax - y || 2 

\\Ax C c-y\\ 

16 

17 

12 

38 

18 

2.8(-12) 

1.9(-13) 

64 

41 

32 

97 

80 

6.6(-12) 

7.7(-12) 

256 

49 

47 

131 

156 

2.8(1 1) 

1 -3(- 1 1) 

1024 

50 

75 

181 

253 

8.7(- 11) 

2.7(-l 1) 

4096 

35 

217 

415 

540 

1.7(-9) 

6.5(1 1 ) 


Table 2: Results for u = 1 . 
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n 


n c 

+ l n c 

«CG 

\\Ax - V lb 

5 

IF 

18 

12 

39 

18 

7-7(- 13) 

5 

64 

51 

24 

93 

80 

9.0(12) 

5 

256 

64 

34 

124 

156 

2.9(11) 

5 

1024 

67 

54 

162 

253 

7.9(11) 

5 

4096 

46 

169 

342 

540 

1.9(-10) 

5 

4096 

109 

73 

237 

540 

4.0(-9) 

10 


Table 3: Results for <7 = 1. Shift jx = r m /6. 


Algorithm 1, nc\ the equivalent number of CG iterations taken by Algorithm 
], rig •+• \nc\ the number of iterations taken by unpreconditioned CG, nca\ 
the residuals achieved by CG and by Algorithm 1. The number of Lanczos 
iterations that were used to estimate eigenvalues varied from 10 to 40. 

In our experience, Tj is extremely accurate, but r m may be several times 
larger than a n . We also found that an underestimate of o„, which results 
in a smaller /x than wc would take given perfect knowledge and hence a bias 
in favor of more iterations at Step 1 and fewer at Step 3, is preferable to an 
overestimate. Thus, we replaced our estimate of a n by 

Oi n ~ 7 m j S 

where 6 > 1 is a parameter of the algorithm. Results are given in Table 3 for 
the same class of matrices as above, with <7=1. Note that with this better 
estimate of the optimal shift, wc are doing far better than with unprecondi- 
tioned CG. 

Strang [16] and Chan [7] have advocated circulant preconditioners. Strang 
takes the circulant C = mo d ,J with c* = , k. = 0, ... ,m where 

m = nf 2. The other diagonals of C are determined by symmetry and the 
requirement that C be circulant. Thus C coincides with A in the central 
half of the diagonals. Chan’s choice is to take C to minimize || A - C|| F , the 
Frobenius norm of the difference, among all circulant C. For this, one takes 
Cfc = ( ka n -k + (n — k)ak)/n. The cost of each iteration is 1.5 times greater 
than the cost of a CG iteration, since solving the circulant preconditioning 
system at each step requires both a forward and an inverse FFT of length n. 


12 




n 

1.5 n s 

1.5 n c 

16 

33 

30 

64 

122 

114 

256 

662 

293 


Table 4: Results for Chan and Strang preconditioners, a — 1. 

We used these techniques for the model problem above, with o — 1. J he 
results, in Table 4, show that neither technique is competitive with ours 
for this problem. For other problems, in particular when a k decays with 
increasing k, we have found them to be superior, as has been predicted by 
Strang, Chan, and Edelman ([17]). 


5 Conclusions 

We have demonstrated a preconditioner for Toeplitz systems that achieves a 
significant speedup when compared with unpreconditioned conjugate gradi- 
ent method and, for certain problems, is considerably better than other previ- 
ously proposed preconditioners. The asymptotic complexity is 0(n logn/c ^ ) 
where k is the condition number of the problem the exponent 1/4 may be 
made smaller at the expense of an increase in the constant factor by applying 
our technique recursively. Compared with the complexity (Bn log 2 n) of other 
superfast methods, we can see that for large, well-conditioned problems the 
present technique may be quite useful. 
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